Radial Spreading of Drift Wave-Zonal Flow Turbulence via Soliton Formation 
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The self- consistent spatiotemporal evolution of drift wave (DW) radial envelope and zonal flow 
(ZF) amplitude is investigated in a slab model |l|. Stationary solution of the coupled partial 
differential equations in a simple limit yields formation of DW-ZF soliton structures, which propagate 
at group velocity depending on the envelope peak amplitude. Additional interesting physics, e.g. 
generation, destruction, collision and reflection of solitons, as well as turbulence bursting can also 
be observed due to effects of linear growth/damping, dissipation, equilibrium nonuniformities and 
soliton dynamics. The propagation of soliton causes significant radial spreading of DW turbulence 
and therefore can affect transport scaling with system size by broadening of the turbulent region. 
Correspondence of the present analysis with the description of DW-ZF interactions in toroidal 
geometry 0, 0] is also elucidated. 



Explaining the size scaling of confinement properties 
in magnetized plasmas is one of the crucial and chal- 
lenging problems of fusion energy research. It has been 
pointed out that turbulence spreading is responsible for 
local turbulence intensity dependence on global equilib- 
rium properties [4|, i.e. the system size, and, thus, for 
the size scaling of turbulent transport coefficients. There- 
fore, the nonlocal character of turbulent intensity plays 
a crucial role in the breakdown of gyro-Bohm scaling of 
turbulent transport and transition to Bohm scaling, as 
observed in several numerical simulations [2, S @] • 

The radial propagation of drift- wave (DW) turbulence 
in tokamak plasmas was first investigated by Garbet Q , 
in the absence of zonal flow (ZF). Turbulence spreading 
was investigated also in Refs. [1, 0]. Later on, using a 
single model equation for the local turbulence intensity, 
Hahm et. al. [10] considered the "minimal problem" for 
turbulence spreading, which is about spatiotemporal dif- 
fusive propagation of a patch of turbulence as a fluctua- 
tion front from an unstable to a stable or a weaker drive 
region. A mean field theory of turbulent transport has 
been developed and extensively studied. By performing 
a Fokker-Planck analysis on the evolution of turbulence 
energy density, or applying quasilinear theory to the wave 
kinetic equation, one can derive a simple equation for the 
mean turbulence energy density. This approach leads to 
a reaction diffusion equation similar to the well-known 
Fisher equation Jj], [l2[ . Giircan et al. [l3[ obtained an 
exact solution for this model, which describes a ballis- 
tic front propagation with speed given by the geometric 
mean of diffusion coefficient and linear growth. In this 
work it was pointed out that ballistic spreading is possible 
even without toroidal coupling effects. A more system- 
atic approach [3] was proposed to explain turbulence 
spreading in terms of nonlinear mode couplings using a 
two field Hasegawa-Wakatani model (kinetic and inter- 
nal energy) recovering the previous one-field model 



in the proper limit, where the fluxes due to nonlinear in- 
teraction are written in the Fick's law form. Analyses of 
turbulence spreading based on solutions of a bi-variate 
Burgers equation [15[ for the evolution of the DW plas- 
mon density were reported in Ref. [16]. Garbet et al. 
also developed a two-field critical gradient model that 
couples a heat equation to an evolution equation for the 
turbulence intensity. It is shown that this model exhibits 
the dual character of turbulent dynamics, diffusive or bal- 
listic, depending on parameters such as the heat flux and 
the wave number. 

In spite of great efforts, the fundamental dynamics of 
turbulence spreading is still not well understood. Al- 
though turbulence is truly a microscopic phenomenon, 
spreading or propagation of turbulence is usually re- 
lated to mesoscale dynamics, e.g. intermittency, for- 
mation of avalanches, transport barriers and other co- 
herent structures, which cannot be described by linear 
excitation and nonlinear wave-wave couplings via triad 
interaction processes only. ZFs are frequently assumed 
to be less or not important at all in the spreading pro- 
cess EH, [HI , based on the argument that large scale 
radially extended eddies are most effective at spreading 
turbulence, while ZFs inhibit spreading by destroying 
these structures H, [Hj. However, slower DW turbu- 
lence spreading, observed in global gyrokinetic simula- 
tions when ZFs are included, has been attributed to the 
suppression of DW intensity by the ZFs [20, |2l| and not 
to their dynamic role. 

In the present work, we study the nonlinear DW-ZF 
interplay in a simple slab geometry in order to elucidate 
the underlying physics mechanisms responsible for tur- 
bulence spreading. A general two-field DW-ZF model is 
derived for the spatiotemporal evolution of the DW ra- 
dial envelope and ZF amplitude, which reduces to previ- 
ous descriptions [2, y, [22( when ZF induced modulations 
on a given DW pump are considered with its sidebands 
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FIG. 1: Collision between two different DW-ZF solitons, 
where A dl = 0.002, k x p s = 0.3 and A d2 = 0.0003, k x p s = 0.3. 

(4- wave). Since the total energy cascades into shorter ra- 
dial wavelengths via the coherent nonlinear DW-ZF mod- 
ulation interaction, the local DW envelope nonlinearly 
steepens and the DW linear dispersion becomes stronger. 
Time scales for nonlinear interaction and linear disper- 
sion eventually become comparable, showing analogies to 
the Langmuir soliton problem. Coherent structures are, 
thus, expected to form, such as DW-ZF solitons which 
will propagate radially. Turbulence spreading may then 
occur via DW-ZF soliton propagation with x ~ t, which 
is faster than any diffusive process. 

The coherent 4- wave DW-ZF modulational interaction 
model [22|] has been applied to study turbulence spread- 
ing in toroidal plasmas, demonstrating that coupled non- 
linear evolution equations for DW radial envelope and 
ZF structures can generally be derived from first princi- 
ples In this work, the same approach is applied in 
a simplified slab geometry [1], where x,y,z corresponds 
to toroidal coordinates r, £, respectively, and the radial 
wave number k x is equivalent to k r = nOkdq/dr [1,0, Q 
of the DW radial envelope. This simplified approach 
also helps elucidating "the subtle differences between the 
slab and toroidal geometries" [l| and yields nonlinear 
equations that are derived from first principles as those 
of [3, Q . This important point allows us to support the 
generality of our results reported hereafter, which do not 
rely on any ad-hoc model assumption for the description 
of nonliner DW-ZF interplay. 

We start from the slab analysis of the electrostatic 
DW-ZF interaction model proposed in Similar to 
the Hasegawa- Mima's model, using two-fluid description 
and quasi-neutrality condition, one can straightforwardly 
derive the DW evolution equation in the form [l|: 

(1 - p 2 s V 2 )d t ^ d - (c 2 s /ni)Vcf) d x z • Vlnno 

-{c 2 s /Vti)V(j> z x z • V0 d + (c 2 a p 2 JSli)V • [V<^ x z 
■VV^ d + V(/) d x z • VV<^] = ; (1) 



where c s = y/T e /rrii, p s = c s /Qi is the ion Larmor radius 
at the sound speed and the scalar potential is normalized 
to T e /e; the ZF potential <p z = ((/>), where (• • • ) repre- 
sents flux surface averaging ((y,z) plane). The last two 
terms on the right-hand side correspond to higher order 
Reynolds stress corrections due to nonlinear polarization 
drift, 0(/c 4 p 4 ), which can be ignored when \kp s \ <C 1. 
Meanwhile, ZF has kg = fcy = 0; thus, electrons do not 
behave adiabatically in the ZF potential. We can de- 
scribe ZFs by the condition of no net radial flux: 

d t V 2 z - (^M)(V • [V0 d x z • VV0 d ]> = . (2) 

As in Refs. [H, 0, Q, we consider a coherent drift wave 
with single toroidal number n , or constant k y in slab ge- 
ometry. Thus, the 2-field coupled set of DW-ZF evolution 
equations are readily cast in the form 

(1 + k 2 - dl)d t (t>d + iu*(x)<j>d = -iC(j> d d x <j> z , (3) 

d t (j) z = iC^ddx^d ~ c ' c -) 5 ( 4 ) 

where uo* — —k y (c 2 s /Vti)dhin/dx is the diamagnetic drift 
frequency, C = k y VLi/uo* is a constant, while space and 
time have been normalized to p s and c^" 1 respectively. 
Note the structural similarity of Eqs.(|3|)-(|1|) with Eqs. (4) 
of [2j] in toroidal geometry. Numerical simulations of the 
above coupled system, given <pd as a Gaussian function 
of x at t = 0, show that DW-ZF can form solitary struc- 
tures, which coherently propagate with given group ve- 
locity (Fig.l). These coherent structures are envelope 
solitons with wavelength of the carrier wave comparable 
to the envelope width, suggesting that turbulence spread- 
ing can be caused by soliton formation due to balance 
between DW dispersion and trapping by nonlinearly gen- 
erated ZFs. 

For the sake of simplicity, we initially ignore linear 
growth/damping and dissipation of both DW and ZF. 
For now, we also take ou* constant. The w*(x) profile in- 
troduces extra effects of finite system size, which will be 
discussed elsewhere. Furthermore, we assume a coherent 
DW form (/>d(x,i) = AdUd(x,t)exp(ik x x — icji), in which 
Ad is the maximum perturbation amplitude, usually 
« 10 -4 — 10 -2 , the normalized envelope function Ud(x : t) 
is chosen to be real and long-scale \d^Ud\ <C ky\ud\, the 
phase cp = k x x — uot describes fast oscillations in time but 
not necessarily in space, k x is the radial wave number 
and cj is the DW frequency. Once the given DW form 
is substituted into Eqs.(|3])-(j4j), the coupled PDEs can be 
rewritten in the form of a nonlinear Schrodinger equation 

(1 + kl)(d t + v g d x )u d + (iuj - dt)d 2 x u d 

-iuu 2 \u d = -iCd x (j) z u d , (5) 
d t (t> z = 2Ck x A 2 u 2 ; (6) 

where v g = —2k x v/(l + fcjj, A = (1 + k\) — and 
k\ = k 2 + k 2 . For constructing a stationary solution, 
we introduce £ = (x — v g t)S and then assume u d (x : t) = 
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FIG. 2: The relation between v g and k x for Ad = 0.003. 



Ud(£) and <p z (x,t) = </>z(£)> sucn that ®t = ~ v g^9^ and 
= 5d%. Here, the small parameter S corresponds to 
the slowly varying envelope scale. Finally, substituting 
<j) z from Eq. (j6]) into Eq. ([5]) , we derive one single ordinary 
differential equation (ODE) for the DW perturbation: 

5 2 u d - \u d + (C 2 /co 2 )(l + k 2 ± )A 2 u d = ; (7) 

where terms ex d t + v g d x cancel by construction and 
d 2 d t Ud is ignored assuming that the envelope transient 
time is much longer than the DW oscillation period, e.g. 
v g S <C uj, which can be justified a posteriori. The above 
second order ODE clearly indicates the competition be- 
tween linear dispersion and nonlinear self-trapping pro- 
cess. When the DW amplitude Ad increases, its envelope 
becomes nonlinearly steeper, i.e. 5 increases; meanwhile, 
the DW dispersion also becomes stronger and tends to 
inhibit the focusing process. Formally, this corresponds 
to equating the three coefficients of v!' d , Ud and u d , i.e. 

S 2 = 1 + k 2 ± - 1/u) = C 2 /(2u 2 )(l + k 2 ± )A 2 (8) 

The DW wave-packet frequency uj is then readily ob- 
tained from the above quadratic equation, i.e. uj = 
[1 + y/l + 2(l + k 2 ± ) 2 C 2 A 2 d ]/[2(l + fcjj] . Note that the 
right-hand side contains the nonlinear frequency shift due 
to finite DW turbulence amplitude. Similarly, the param- 
eter S can also be determined as 
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kl)/2C(A d /uo) 



(9) 



This derivation is subject to our a priori assumption 
S <C uVg 1 , which guarantees that the DW oscillation uj 
occurs on the fastest timescale. Substituting 5, uo,v g as 
functions of k x and Ad, this assumption is equivalent to 
k x kyVLi/uo*A d < + k 2 _)/2uj(k x ,A d ). Finally, given 
Eq.®, Eq.([7j) becomes a dimensionless ODE governing 
the stationary envelope function Ud, 



This ODE is analogous to that of an oscillator in the 
so called "Sagdeev potential" $>(ud) = (—u 2 + 
whose solution can be written as hyperbolic secant 
function Ud(£) = Sech(^), when appropriate bound- 
ary conditions are imposed, viz. — > at |<^"| — > 
oo. Meanwhile, the ZF solution is obtained straightfor- 
wardly by integrating Eq.(|6|) once, such that 4> z (£) = 
f2k x C/(v g 5)A*Sech 2 (€)d£ = ^^2(1 + k 2 ± ) A d Tanh(0 
which statisfies the causality constraint, i.e. d^<p z —> 
when |^| oo for any initially localized DW turbulence. 
The expressions for DW and ZF in the laboratory frame 
are 

ct> d {x,t) = A d Sech[5(x + r ?^t)]e ik * x - i " t , (11) 



2k x uo 



v!' d -u d + 2u 6 d 



0. 



(10) 



<l> z (x,t) = v /2(l + fc2)A d Tanh[5(a;+ j^r*)] • (12) 

From Eqs.jni) and (O, it generally follows that ZF 
potentials have radially moving structures of hyperbolic 
tangent shape; meanwhile, E z = —d(j) z /d x manifests it- 
self as scalar-potential wells in the background plasma 
and trap the corresponding DW packets. Figure [T] shows 
the spatiotemporal evolution of two counter-propagating 
DW-ZF solitons, which are solutions of the original cou- 
pled PDEs, given k y = 0.3, Qi/u = 100. For consistency 
with our analytic approach, we have chosen initial k x and 
Ad to satisfy the a priori assumption S <C wVg 1 . Fur- 
thermore, we assumed no oo* equilibrium variation and 
no growth/damping and dissipation. Note that the two 
envelope solitons remain unchanged in both real and k 
space after the collision, although the dynamics during 
the collision can be quite complicated. This is one of the 
soliton essential features. 

The radial propagation velocity of DW-ZF solitary 
structures, v g , depends on both the radial wave number 
and the DW amplitude. It is different from the linear 
group velocity, v' g = duoi/dk Xl which is determined by 
k x through the linear dispersion relation only. There- 
fore the solution of Eqs. ([TTj) - ([T2]) gives a two-parameter, 
k x and Ad, family of solitons. Figure [2] shows the re- 
lation between v g and k x for small initial amplitude 
Ad = 0.003. Numerical and analytical results agree well 
when k x < 0.4. The discrepancy when k x > 0.4 origi- 
nates from the breaking of the a priori assumption that 
d 2 dt4>d can be ignored in Eq.([5]). Moreover, when Ad in- 
creases to about 10 -2 , the analytical result is no longer 
valid either, since the ignored term 0(v g 5 3 ) modifies the 
solution at larger k x or Ad, according to Eq.©. 

Our numerical simulation results for Ad > 0.01 show 
that the dominant asymptotic (t —> oo) DW turbulence 
behavior is still of soliton type and the propagation veloc- 
ity v g increases with the DW amplitude Ad- We observe 
that the DW radial wave number no longer corresponds 
to its initial value but rather to 5, which is mainly deter- 
mined by the amplitude Ad alone. There seems to be a 
transition from a 2-parameter to a 1-parameter family of 



4 



(a) Radial Profile 
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(b) DW Turbulence Amplitude 
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FIG. 3: DW turbulence spreading in the presence of 
growth/damping, dissipation and finite system size effects. 



soliton solutions of the coupled system. If the amplitude 
becomes even larger, e.g. ^ 0.02, the initially local- 
ized DW-ZF soliton breaks into many pieces in the form 
of solitons and wave-trains; similar to Gardner's work [23[ 
on Korteweg-deVries equation, in which it is shown that a 
localized but otherwise arbitrary initial perturbation will 
generate a conventional wave- train, quickly destroyed by 
dispersion, and a finite number of solitons, which char- 
acterize the asymptotic solution. 

We studied the DW-ZF initial value problem in more 
general cases as well, i.e. in the presence of linear 
growth/damping, dissipations, and variation of equilib- 
rium profiles. Figure [3] shows the evolution of DW turbu- 
lence out of initial random noise, with strong DW growth 
rate 7d(0) = 0.1, uniform ZF damping rate j z = 0.075 
and L p = 150, which represents the system size. Dissi- 
pations are also included. The drift frequency uj*(x) has 
Gaussian shape centered at x = (Fig. [31(a)); DW tur- 
bulence is linearly unstable in the central region (\x\ < 
S0p s ) but is damped in the outer region (\x\ > 80p s ), 
while the ZF is uniformly damped. Figure [31(b) clearly 
shows formation and propagation of solitons, which how- 
ever exhibit more complicate dynamic behaviors; for in- 
stance, growing amplitudes, slowing down of propaga- 



tion speed, soliton breaking, turbulence bursting and 
more. Since coupled PDEs generally describe infinite- 
dimensional dynamical systems, DW turbulence dynam- 
ics appears mostly chaotic in the corresponding param- 
eter space, (7d,7z). Solitons may bounce back at their 
turning points, possibly enhancing nonlinear interactions 
inside the turbulent region and impacting the size scaling 
of turbulence transport. Figure [31(b) also demonstrates 
that the nonlinearly saturated turbulence has spread into 
a much broader region than that of its linear mode struc- 
ture. 

In summary, we have demonstrated the novel result 
that coherent structures such as solitons can be con- 
structed self-consistently in a two-field DW-ZF model 
and cause significant radial turbulence spreading in a slab 
plasma. We have also shown the structural analogy of the 
underlying coupled PDEs for the nonlinear evolution of 
DW radial envelope and ZF amplitude with the corre- 
sponding equations derived in toroidal geometry [!, 0], 
demonstrating the generality of the present results and 
the possibility of readily extending them in future works. 
The size scaling of DW turbulence will also be discussed 
in detail in a separate work. 
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